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^  ABSTRACT 

An  implicit  finite  difference  model  for  predicting  flood  routing  is  applied  to 
the  lower  Columbia  River,  where  tidal  forcing  causes  flow  reversals  interacting 
with  upstream  dam  flow  during  small  river  flow  periods.  The  model  is  one¬ 
dimensional,  unsteady,  including  lateral  inflow  and  variable  bed  friction  for 
different  channel  sections.  A  comparison  of  stages  at  six  stations  was  made  for  a 
sensitivity  analysis.  The  analysis  used  a  total  of  2209  hours  of  simulated  river 
stages.  /-  ^  -T) 

Downstream  boundary  changes  of(^b.5  feet  andd&.O  feet  were  made  to  the 
Astoria  tide  stages.  Model  simulations  showed  that  70%  of  the  tide  difference 
appears  at  Vancouver  and  Portland,  80%  at  St.  Helens,  85%  at  Longview,  93%  at 
Wauna  and  95%  at  Skamokawa.  Varying  the  upstream  boundary  condition 
(Bonneville  Dam  discharges)  by  <$10%  and  were  markedly  different  from 
the  downstream  boundary  changes.  Upstream,  where  the  tide  influence  is  weakest, 
the  tidal  cycle  is  more  likely  to  be  “washed  out”  by  the  higher  flows  of  the 
Columbia.  Also  these  changes  fluctuated  with  the  tide  cycle.  Downstream  stations 
did  not  show  such  differences  because  of  the  larger  cross  section  areas  of  the 
Columbia  River  nearer  the  mouth  and  the  proximity  to  the  downstream  boundary 
condition. 

The  river  system  was  calibrated  in  a  downstream  to  upstream  direction  and 
used  a  total  of  606  hours  of  observed  river  stages.  Three  periods  with  distinct  river 
flow  conditions  were  used  in  the  calibration.  Regression  analyses  of  the  computed 
residual  values  for  each  of  the  stations  gave  correlation  coefficients  (r2)  less  than 
0.360.  However,  cross  correlations  between  residual  and  computed  stages  showed 
that  the  two  series  were  highly  sinusoidally  correlated  for  all  stations.  A  spectral 
estimation  of  the  residuals  exhibited  very  stroiig  peaks  at  frequencies  of  0.081  hr1 
(12.3  hrs),  0.042  hr1  (24.0  hrs)  and  subsequent  harmonics  of  these  frequencies. 
The  residual  components  are  strongly  associated  with  the  tidal  cycle. 
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I.  INTRODUCTION 


A.  BACKGROUND  INFORMATION 

The  accuracy  of  forecasted  stages  is  critical  for  both  navigational  safety  and 
commerce  on  any  river  system.  This  is  especially  true  on  the  lower  Columbia 
River  as  ship  traffic  increases.  The  Northwest  River  Forecast  Center  (NWRFC, 
part  of  the  National  Oceanic  and  Atmospheric  Administration,  NOAA)  presently 
uses  the  Dynamic  Wave  Operational  Model  (DWOPER)  to  forecast  hourly  river 
stages  on  the  lower  Columbia  River  from  the  Vancouver- Portland  area  to  Astoria, 
Oregon  (Figure  1).  This  model  was  developed  by  Fread  (1976),  of  the  National 
Weather  Service  Office  of  Hydrology. 


Figure  1.  DWOPER  Simulation  Area  (not  to  scale) 
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Currently,  DWOPER  is  being  used  on  the  Lower  Mississippi  River  System 
(470  km  or  292  miles),  the  Mississippi-Ohio-Cumberland-Tennessee  River  System 
(633  km  or  393  miles)  as  well  as  the  Columbia-Willamette  River  System  (230  km 
or  143  miles).  It  is  the  Columbia-Willamette  System  that  has  the  most  pronounced 
tidal  influences. 

The  Columbia-Willamette  River  system  is  a  major  waterway  in  the  northwest 
United  States  and  has  a  shallow  slope  (0.01 1  m/km)  below  Bonneville  Dam.  Tidal 
effects  on  the  river  extend  as  far  upstream  as  the  Bonneville  Dam  tailwaters,  230 
km  (143  mi)  from  the  ocean,  during  periods  of  low  flow.  Flow  reversals  can  occur 
as  far  upstream  as  the  Vancouver-Portlar.d  area,  170  km  (106  mi)  from  the  ocean 
(Fread,  1976,  1978).  As  an  additional  part  of  the  river  optimization  process  for 
ships,  the  U.S.  Army  Corps  of  Engineers  has  deepened  the  bar  crossing  at  the 
Pacific  Ocean  to  16.8  m  (55  ft)  and  maintains  a  12.2  m  (40  ft)  channel  to  the  Van- 
couver-Portland  port  area,  a  distance  of  193  km  (120  miles)  from  the  ocean  to  ac¬ 
commodate  larger  draft  vessels  (Rooks,  1986). 

B.  OBJECTIVE 

This  thesis  examines  the  DWOPER  computer  model  (program  version 
4/11/79)  to  find  parameter  values  needed  to  best  match  computed  stages  with 
observed  stages.  Calibration  periods  come  from  the  1987  Water  Year  (October 
1986  to  September  1987).  These  comparisons  are  made  at  six  stations  for  three 
different  time  periods. 

Confidence  intervals  for  forecasted  stages  are  developed.  Residual  errors  be¬ 
tween  model  output  and  observed  stages  are  examined  as  well.  The  modelling  of 
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these  residual  errors  and  applying  confidence  intervals  will  aid  in  more  accurate 
river  stage  forecasts  and  improve  the  DWOPER  application. 

C.  AREA  DESCRIPTION 

The  Columbia  River  dominates  the  Pacific  northwest  from  its  source, 
Columbia  Lake  in  Canada's  Selkirk  Mountain  Range  at  809  m  (2650  ft)  above  mean 
sea  level.  It  is  the  fourth  largest  river  in  North  America  with  a  course  length  of 
1,953  km  (1,214  mi)  and  a  drainage  area  of  671,000  km2  (259,000  mi2).  This  river 
system  is  important  to  agriculture,  fisheries,  power  production,  recreation  and 
commerce  in  the  surrounding  area.  The  Columbia  is  one  of  the  world's  most 
managed  rivers.  Portland,  Oregon,  at  the  confluence  of  the  Columbia  and 
Willamette  rivers,  is  the  third  busiest  port  on  the  west  coast  of  the  U.S.  in  terms  of 
total  tonnage  shipped  per  year.  There  are  38  ports  that  are  served  by  the  Columbia 
system  (including  the  Snake  River)  with  the  farthest  inland  being  Lewiston,  Idaho. 

The  Columbia  River  Management  Report  for  1987  stated  that  precipitation  and 
runoff  was  significantly  below  normal  in  1987  for  the  Columbia  River  Basin. 
Winter  precipitation  was  generally  60%  to  80%  of  normal.  The  spring  runoff  of 
the  Columbia  at  The  Dalles,  Oregon,  was  about  70%  of  normal.  The  regulated 
peak  flow  was  7,950  m3/sec  (284,000  ft  3/sec)  compared  to  a  bankful  flow  of 
12,600  m3/sec  (450,000  cfs).  This  was  the  lowest  runoff  since  1977  and  it  ranks  as 
the  fifth  lowest  compared  to  the  1928-1978  period  of  record.  Snake  River  runoff 
was  only  46%  of  normal  and  ranked  the  second  lowest  during  this  same  50  year 
period  (Columbia  River  Management,  1987  ). 

Despite  the  low  runoff,  reservoirs  in  the  power  system  refilled  due  in  part  to 
the  Bonneville  Power  Administration  (BPA)  conservative  marketing  policy. 
However,  irrigation  reservoirs  were  at  very  low  levels  and  had  little  carryover  for 
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1988.  The  dry  weather  in  the  Cascade  region  caused  several  reservoirs  in  the 
Willamette  basin  not  to  refill  in  the  spring.  August  and  September  were  unusually 
dry  in  some  areas  and  resulted  in  greater  than  usual  reservoir  drawdown  by  the  end 
of  the  1987  water  year  (Columbia  River  Management  Group,  1987  ). 

The  Pacific  Ocean  tides  are  characterized  as  semi-diurnal  (two  highs  and  two 
lows  per  day)  with  a  mean  range  of  2.4  to  3.0  m  (8  to  10  ft).  No  unusual  tidal  phe¬ 
nomena  (i.e.  tsunami)  occurred.  The  Tongue  Point  gauge  at  Astoria,  Oregon  is  a 
primary  station  for  the  National  Ocean  Service  (NOS)  with  more  than  19  years  of 
continuous  record  and  is  a  published  prediction  site.  Values  for  this  study  used 
observed  Tongue  Point  records.  These  values  were  obtained  from  the  Tidal 
Analysis  office  of  National  Ocean  Service  (NOS  Pacific  Tide  Tables,  1987). 

D.  PRESENT  USES  AND  USERS 

The  most  dramatic  use  of  DWOPER  on  the  lower  Columbia  River  occurred 
during  the  May,  1980,  Mount  Saint  Helens  volcano  eruption  in  Washington  state. 
The  NWRFC  used  DWOPER  to  aid  the  U.S.  Coast  Guard  in  timing  vessel  traffic 
arrivals  at  the  Longview  constriction  for  safer  passage  (Orwig,1980).  Vessels 
arrived  at  or  near  high  water  to  have  maximum  water  depth  over  the  constriction 
created  by  the  mud  and  debris  flows  from  the  volcano  into  the  Columbia  River.  The 
use  of  the  model  during  this  crisis  showed  the  benefits  of  forecasting  river  stages. 

In  February  1984,  a  Sea  Use  Marine  Services  conference  was  held  in  Portland, 
Oregon.  Representatives  of  the  Port  of  Portland  (POP),  U.S.  Army  Corps  of  En¬ 
gineers  (COE),  National  Weather  Service  (NWS)  and  National  Ocean  Service 
(NOS),  both  part  of  the  NOAA,  as  well  as  public  and  private  organizations  at¬ 
tended.  Two  action  items  resulting  from  this  meeting  were  the  need  to  establish  an 
automatic  real-time  river  gauging  network  and  the  need  to  develop  an  hourly  fore- 
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cast  model  for  the  tidal  stages  on  the  lower  Columbia  River.  The  real-time  data  and 
forecasts  would  be  extremely  valuable  to  the  commercial  river  traffic  and  port  au¬ 
thorities  along  the  lower  Columbia.  It  was  at  this  time  that  preliminary  work  for 
using  DWOPER  for  forecasting  these  hourly  stages  on  the  river  began.  Eventually 
it  was  decided  that  the  Port  of  Portland  would  operate  and  maintain  the  gauging 
stations  and  the  automatic  data  collection  system.  The  NWRFC  agreed  to  prepare 
and  issue  a  three  day  forecast  for  the  six  sites  on  the  river  (Orwig  et  al.,  1986). 

The  Port  of  Portland  and  the  Columbia  River  Pilots  Association,  through  the 
local  merchants  exchange  group,  have  been  distributing  the  three  day  forecasts 
from  the  NWRFC  to  area  shipping  interests.  David  Nesset,  director  of  marine  ser¬ 
vices  for  the  Port  of  Portland  has  stated,  “A  foot  of  difference  (0.3048  m)  on  a  ship 
is  a  whole  bunch  of  cargo  and  a  whole  bunch  of  revenue.”  He  estimated  the  extra 
revenue  at  $10,000  to  $20,000  for  most  grain  ships.  The  extra  feet  are  helping 
change  the  image  of  the  Columbia  River,  which  used  to  be  thought  capable  of  han¬ 
dling  only  smaller  draft  ships.  Lawrence  H.  Bogle,  plant  manager  of  the  Peavey 
Grain  Company  of  Kalama,  Washington  explained  that  if  the  Columbia  can 
accommodate  the  so-called  Panamax  ships,  the  maximum  size  of  ship  that  can 
transit  through  the  Panama  Canal,  it  would  mean  that  the  local  river  shipping 
system  “can  compete  with  anybody,  including  the  Mississippi  River.”  Panamax 
ships  have  drafts  up  to  13.1  m  (43.0  ft).  In  1984, 28%  of  the  cargo  ships  had  a  draft 
of  more  than  39  feet  and  the  percentage  of  ships  with  larger  drafts  increased  to  68% 
in  1985  (Rooks,  1986). 
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IL  DWOPER  MODEL  DESCRIPTION 


A.  BACKGROUND  INFORMATION 

Unsteady  flow  in  rivers,  reservoirs  and  estuaries  is  caused  by  the  motion  of 
long  waves  from  floods,  tides,  storm  surges  and  reservoir  releases.  Since  this  flow 
can  be  considered  as  one-dimensional,  the  acceleration  and  velocity  components  of 
the  flow  in  the  transverse  and  vertical  directions  can  be  neglected,  being  small 
compared  to  the  acceleration  and  velocity  components  in  the  direction  of  river 
flow.  This  allows  the  flow  motion  to  be  described  by  the  one  dimensional  Saint- 
Venant  equations  of  mass  and  momentum  conservation  (Fread,  1976). 

The  Dynamic  Wave  Model  (DWOPER)  allows  for  dynamic  routing  of  rivers 
in  such  cases  as 

•  upstream  movement  of  tides  and  storm  surges, 

•  reservoir  backwater  effects, 

•  flood  waves  on  shallow  sloped  rivers,  or 

•  abrupt  waves  from  reservoir  releases  or  dam  failures. 

B.  THEORETICAL  DEVELOPMENT 

Fread  (1976)  developed  the  Dynamic  Wave  model,  assuming  the  following 
conditions: 

•  depth,  velocity  and  acceleration  vary  only  in  the  longitudinal  axis  of  the  river, 

•  velocity  is  constant  for  any  cross-sectional  area  of  the  river  and  that  the  water 
surface  is  horizontal  to  any  section  perpendicular  to  the  longitudinal  axis  of  the 
river, 

•  the  longitudinal  axis  of  the  river  is  approximately  a  straight  line, 

•  channel  bottom  slope  is  small, 
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•  no  erosion  or  deposition  along  the  channel  bottom  occurs, 

•  flow  is  gradually  varied  with  the  hydrostatic  pressure, 

•  the  Manning  equation  describes  the  resistance  effects  and 

•  the  fluid  is  incompressible  and  homogeneous. 

The  model  equations  are: 


conservation  of  mass: 

30  cKA+Aq) 

(1) 

(Note:  The  conservation  of  mass  for  steady  state  flow  is  expressed  as  Q  = 

VA). 

conservation  of  momentum: 

(2) 

friction  slope: 

-2  IQI  Q 
*  2,2  A2  R4/' 

(3) 

Ke  d(%2 

energy  loss  slope: 

(4) 

where, 

Q  =  discharge,  A  =  cross  sectional  area,  x  =  distance  along  the  channel, 

Aq  =  off  channel  cross  sectional  area  (flow  velocity  is  considered  negligible), 

q  =  lateral  inflow  or  outflow,  t  =  time, 

vx  =  velocity  of  lateral  inflow  in  the  x  direction, 

g  =  gravity  acceleration  constant,  B  =  channel  top  width,  Sf  =  friction  slope, 
K«  =  contraction/expansion  factor. 
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The  Ke  term  is  used  to  account  for  local,  small  scale  head  losses  that  may  form 
in  the  river  flow  at  abrupt  cross  section  changes  along  the  river  such  as  a  bridge 
support  column.  For  a  surface  cross  section  of  width  b,  the  affect  is  usually  limited 
to  a  distance  of  about  10-b  (personal  conversation,  Fread).  The  Ke  term  was  set  to 
zero  for  this  study. 

Equations  (1)  and  (2)  are  nonlinear  partial  differential  equations  that  can  be 
solved  by  explicit  or  implicit  finite  difference  techniques.  Explicit  techniques  are 
generally  not  well  suited  for  application  to  long  period,  unsteady  flow  situations 
such  as  flood  or  tide  waves.  These  techniques  are  restricted  by  numerical 
instability  reasons  to  very  small  computational  time  steps  (the  order  of  a  few 
minutes).  This  would  be  a  poor  and  inefficient  use  of  computer  resources.  On  the 
other  hand,  implicit  finite  difference  techniques  have  less  restriction  on  time  step 
size,  except  for  accuracy  considerations. 

The  “weighted  four  point”  method  used  in  DWOPER  was  first  used  by 
Preissmann  (1961).  This  method  has  the  advantage  of  using  unequal  distance  steps 
and  controlled  stability-convergence  properties  (Fread,  1978).  A  continuous  x,t 
region  represented  by  a  grid  of  discrete  points  allows  for  solutions  of  both  h  and  Q 
at  equal  or  unequal  intervals  of  Ax  and  At  along  the  x  and  t  axes  where  i  is  the  x 
position  and  j  the  time  position  (see  Figure  2).  The  time  derivative  for  any 
variable,  k ,  is  approximated  by  equation  (6). 


3k 

3T 


K1 


+  k 


j+1  . 

i+1 


kJ.  - 

l 


J 
i+1 
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2At 


(6) 
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The  spatial  derivative  is  approximated  by  a  finite  difference  quotient 
positioned  between  two  adjacent  lines  (see  Figure  3)  according  to  the  weighting 
factors  0  and  (1*8),  such  that: 


dk 

ST 


8 


Ax 


+ d-e) 


Ax 


(7) 


with  the  non-derivative  terms  approximated  by: 


k«8 


tkji+1  • k  r+‘  ]  „  tk  \  -  *  ii] 


+  d-8) 


(8) 


When  8  equals  1.0,  a  fully  implicit  scheme  is  formed  (Baltzer  and  Lai,  1968) 
and  when  8  equals  0.5,  a  box  scheme  is  formed  (Amein  and  Fang,  1970).  Fread 
(1975)  studied  the  stability  and  convergence  properties  and  concluded  that  the 
accuracy  decreases  as  8  departs  from  0.5  and  approaches  1.0.  This  effect  is  more 
pronounced  as  time  step  size  increases.  The  model  version  used  in  this  study  allows 
6  to  be  an  input  variable  with  0.55  used  as  the  default  value.  This  value  is  used  to 
minimize  loss  of  accuracy. 

Substitution  of  equations  (6),  (7),  (8)  yields  two  algebraic  equations  that  are 
nonlinear  with  respect  to  h  and  Q  at  the  nodal  points  on  the  j+1  time  line.  All  terms 
from  the  j  time  line  are  known  from  either  initial  conditions  or  previous 
calculations.  The  initial  conditions  are  values  of  h  and  Q  at  each  computational 
point  or  node  along  the  x-axis  for  the  j+1  time  line. 

The  two  algebraic  equations  cannot  be  solved  directly  since  there  are  four 
unknowns  (h  and  Q  at  points  i  and  i+1  on  the  j+1  time  line).  If,  however,  similar 
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equations  are  formed  for  each  of  the  N-l  Ax  reaches  (N  equals  the  number  of  cross 
section  computational  points  on  the  river),  a  total  of  (2N-2)  equations  with  2N 
unknowns  results.  These  equations  along  with  the  upstream  and  downstream 
boundary  conditions  yield  a  system  of  2N  nonlinear  equations  with  2N  unknowns, 
which  are  solved  by  the  Newton-Raphson  iteration  procedure  method  (Fread, 
1976).  The  tolerance  of  the  residuals  produced  by  the  iteration  method  can  be 
specified  as  a  model  input  parameter  and  was  0.001  in  this  study. 

The  Dynamic  Wave  model  allows  initial  conditions  to  be  obtained  from: 

•  estimated  stages  and  discharges  at  each  cross  section, 

•  observed  stages  at  each  cross  section  where  a  gauge  may  be  located  with 

1)  stages  at  intermediate  cross  sections  linearly  interpolated  within  the  model, 

2)  observed  discharges  from  the  upstream  boundary,  or 

3)  with  downstream  discharges  found  by  summing  flows  from  the  upstream 
boundaries  of  the  main  river  and  any  dynamic  tributaries  plus  any  lateral 
inflows, 

•  saved  values  of  stages  and  discharges  from  previous  simulations,  and 

•  assumed  steady  state  flow  to  obtain  discharges  and  backwater  computations  to 
obtain  stages  (Fread,  1978). 

In  each  case,  the  unsteady  flow  equations  are  solved  for  several  time  steps  using 
the  initial  conditions  together  with  boundary  conditions  which  are  held  constant 
during  the  time  steps.  This  allows  the  errors  in  the  initial  conditions  to  be  dampen 
out  which  results  in  the  initial  conditions  being  more  nearly  error  free  when  the 
actual  simulation  begins  and  transient  boundary  conditions  are  used. 
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Figure  2.  Numerical  Grid  Configuration  and  Boundary  Conditions 


C.  BOUNDARY  CONDITIONS 


Boundary  conditions  must  be  specified  in  order  to  solve  the  Saint- Venant 
equations.  In  most  unsteady  flow  problems,  the  unsteady  disturbance  is  introduced 
into  the  flow  at  the  boundaries  or  extremities  of  the  river  system.  DWOPER  can 
accommodate  either  known  h  or  Q  [hi(t)  or  Qi(t)]  as  an  upstream  condition.  The 
model  can  also  accommodate  either  known  h  or  Q  as  a  downstream  condition  [hu(t) 
or  QnCO]  or  use  a  rating  curve  option  (flow  versus  stage  relationship).  If  the  rating 
curve  option  is  used,  the  rating  may  be  single-valued  and  read  in  as  tabular  values. 
Intermediate  points  are  linearly  interpolated.  The  rating  may  also  be  a  loop  rating 
curve  generated  internally  from  cross  section  and  roughness  coefficients  for  the 
downstream  extremity  that  have  already  been  entered  into  the  program  and  the 
instantaneous  water  surface  slope  at  the  previous  time  step.  However,  the  rating 
curve  option  was  not  used  in  this  study.  Ocean  tidal  affects  on  the  river  makes  any 
flow-stage  relation  changeable  with  time  (Fread,  1978). 

Boundary  conditions  in  this  study  are  Bonneville  Dam  discharge,  Willamette 
River  flow  below  the  Estacada  River  confluence  and  the  Pacific  Ocean  tides  at 
Tongue  Point,  Oregon  (near  Astoria,  Oregon).  Flow  data  are  instantaneous  values 
at  one  hour  time  intervals.  If  the  Willamette  River  observed  flow  record  had 
irregularly  spaced  data,  intermediate  flow  values  were  linearly  interpolation  from 
the  observed  data.  There  was  no  missing  data  from  Bonneville  Dam  or  the  Tongue 
Point  gauge. 

A  sample  printout  from  the  DWOPER  model  is  shown  in  Figure  4.  Both 
computed  flow  (Q-FCST)  and  computed  stage  (H-FCST)  are  displayed  with  the 
observed  stages  for  Skamokawa,  Washington. 
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RIVER  I,  STATION  21  SKAMOKAWA  FLOOD  STAGE  -  14.00  FEET 
*  COMPUTED  + - OBSERVED  STAGE (FT) 


TIME  -2. 

0. 

2. 

4. 

6. 

8.  Q-FCST 

H-FCST 

H-OBS 

10/16/13. 

$ 

*+.-260.84 

7.67 

7.86 

10/16/14. 

$ 

+  .  -11.91 

7.33 

7.41 

10/16/15. 

$ 

+ 

.  311.14 

6.02 

6.09 

10/16/16. 

$ 

.  493.99 

4.54 

4.63 

10/16/17. 

S 

+ 

.  508.14 

3.23 

3.18 

10/16/18. 

$ 

+ 

.  483.79 

1.94 

1.96 

10/16/19. 

$  *  + 

.  486.42 

0.61 

0.96 

10/16/20. 

*s  + 

.  416.19 

-0.27 

0.32 

10/16/21. 

$  *  + 

.  226.16 

0.44 

0.72 

10/16/22. 

$ 

.  *  + 

.  27.99 

2.18 

2.37 

10/16/23. 

$ 

*  + 

.-144.99 

4.11 

4  .42 

10/17/  0. 

$ 

*  + 

.-236.32 

5.62 

6.02 

10/17/  1. 

$ 

*+  .-257.95 

6.74 

7.01 

10/17/  2. 

5 

*+  .-106.87 

6.97 

7.13 

Figure  4.  Example  Hydrograph  Output  from  DWOPER 


D.  OTHER  MAJOR  FACTORS 

The  model  uses  a  total  of  29  separate  river  reaches  (24  for  the  Columbia  River 
and  5  for  the  Willamette  River),  each  with  cross  sectional  information  and  river 
mile  location  and  Manning  n  (roughness  coefficient)  values.  An  explanation  of 
each  of  these  items  follows. 

1.  Cross  Sectional  Areas 

Cross  sections  can  be  entered  as  either  regular  or  irregular  shaped  areas. 
Each  is  read  in  as  tabular  values  of  channel  width  and  elevation,  together  giving  a 
piece-wise  linear  relationship.  A  low  flow  area  which  can  be  zero  is  used  to 
describe  the  cross-section  below  the  minimum  elevation  inputted.  During  the 
solution  of  the  Saint  Venant  equations,  areas  or  widths  associated  with  a  particular 
water  surface  elevation  are  linearly  interpolated  from  the  piece-wise  linear 
relationships  inputted.  Cross  sections  from  gauging  station  locations  are  always 
used  as  computational  points  in  the  x-t  plane.  The  average  cross  sectional  area  is  a 
weighted  average  of  the  cross  sectional  properties  of  the  intervening  reach.  The 
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distance  between  each  cross  section  is  used  as  the  weighting  factor.  Observed 
gauging  stations  and  cross  section  locations  used  in  this  study  are  in  Table  1 . 

2.  Lateral  Inflow  Sources 

The  dynamic  wave  model  allows  lateral  inflow  tributaries  (non-dynamically 
routed  tributaries)  that  account  for  the  q  term  in  equations  (1)  and  (2).  These  flows 
are  considered  independent  of  the  main  river  flow  and  are  read  in  the  program  at 
constant  or  varying  time  intervals.  Any  location  along  the  main  river  or  a 
dynamically  routed  tributary  may  be  specified.  Flow  units  are  cubic  feet  per 
second  (cfs)  with  outflows  being  negative  values.  Flows  at  other  times  are  linearly 
interpolated  between  the  appropriate  inputted  time  values.  Lateral  inflow  locations 
are  shown  in  Table  2.  The  Vancouver  local  inflow  term  in  Table  2  simply  refers  to 
the  cumulative  inflows  from  all  streams  and  creeks  between  Bonneville  Dam  and 
Vancouver,  Washington 

3.  Friction  terms 

The  Manning  n  roughness  coefficient  is  used  to  describe  the  friction  along 
the  river  bed  caused  by  vegetation,  obstructions,  bend  effects  or  eddy  losses  and  is 
defined  for  each  group  of  reaches  bounded  by  gauging  stations  (see  Figure  5). 
These  roughness  coefficients  can  be  entered  either  as  a  function  of  discharge  or 
stage  height.  Simulations  can  vary  markedly  with  small  changes  in  the  roughness 
coefficients  and  any  trail  and  error  type  calibration  requires  a  knowledge  of  how 
this  input  data  will  affect  model  output. 

In  this  study,  all  roughness  coefficients  are  functions  of  the  river  flow  and 
not  stage.  The  n  values  are  read  into  the  program  as  tabular  values  for  a  specific 
flow  with  intermediate  values  linearly  interpolated.  All  input  parameters  used  in 
this  study  are  listed  in  the  appendix. 
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TABLE  1.  GAUGING  STATION  AND  CROSS  SECTION  LOCATIONS 


Columbia  River 


Willamette  River 


tm 

Mile 

Description 

mm 

Mile 

i 

145.5 

Bonneville  Dam 

i 

24.2 

Gladstone,  OR 

2 

143.3 

Ives  Island 

2 

18.6 

Milwaukie,  OR 

3 

141.1 

Warrendale,  WA 

3 

12.8 

Portland,  OR 

4 

132.0 

Sand  Island 

4 

6.4 

Saint  Johns,  OR 

5 

122.9 

Washougal,  WA 

5 

0.0 

Columbia  Conf. 

6 

114.7 

East  Government  Island 

7 

106.5 

Vancouver,  WA 

8 

103.8 

East  Hayden  Island 

9 

101.6 

Willamette  River* 

10 

101.1 

Willamette  River** 

11 

92.5 

Bachelor  Point 

12 

87.0 

Lewis  River* 

13 

86.9 

Lewis  River** 

14 

86.1 

Saint  Helens,  OR 

15 

75.0 

Kalama,  WA 

16 

67.6 

Cowlitz  River* 

17 

67.5 

Cowlitz  River** 

18 

66.1 

Longview,  WA 

Gauged  stations  are  in  italics 

19 

53.9 

Bradbury  Slough 

20 

41.6 

Wauna,  OR 

*  _ 

Upstream  Confluence 

21 

33.7 

Skamokawa,  WA 

**  = 

Downstream  Confluence 

22 

30.2 

Three  Tree  Point 

23 

23.4 

Miller  Sands 

24 

17.5 

Tongue  Point 

TABLE  2.  LATERAL  INFLOW  LOCATIONS 


Columbia  River 


□ 

Mile 

Description 

u 

jMP£9r£ 

Vancouver  Local 

H 

1 

Lewis  River 

H 

Cowlitz  River 

15 


Boundary 

condition 

Bonneville  Dam 

© 

Boundary 

condition 

Local 

Vancouver 

© 

Lewis  River 

© 

© 
l . 

f 

1  Portland 

St.  Helen* 

Willamette  River 
below 

© 

Estacada  River  confluence 

Cowlitz  River 

Longview 

© 

Wauna 

© 

Boundary 

condition 

© 

Skamokawa 

City 

Astoria  ” 

=  Reach  Number 

=  Observed  Stage  Location 

=  Lateral  Inflow 

PACIFIC  OCEAN 

Figure  5.  Schematic  Diagram  of  the  Lower  Columbia  River 
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III.  MODEL  ANALYSIS 


A.  COMPARISON  DATA 

As  with  any  mathematical  model,  the  accuracy  of  the  model  results  can  only  be 
as  good  as  the  information  that  goes  into  the  program.  In  this  study  of  DWOPER, 
data  for  all  lateral  flows  and  boundary  conditions,  except  Astoria  tides,  were 
obtained  from  the  Northwest  River  Forecast  Center,  Portland,  Oregon.  Ocean  tide 
values  were  obtained  from  the  NOAA’s  National  Ocean  Service  (NOS)  tidal 
analysis  office.  Data  obtained  from  NOS  must  be  adjusted  to  the  datum  the  user  is 
interested  in  using.  Figure  6  shows  different  gauge  readings  for  the  same  water 
level  using  three  different  datums;  Mean  Lower  Low  Water  (MLLW),  Mean 
Higher  High  Water  (MHHW)  and  the  National  Geodetic  Vertical  Datum  (NGVD) 
of  1929.  In  this  study,  all  gauges  are  referenced  to  NGVD  by  means  of  the 
Columbia  River  Datum. 

The  Columbia  River  Datum  is  a  sloping  line  that  connects  the  zero  gauge 
readings  of  the  observation  stations  (see  Figure  7)  and  represents  the  approximate 
lowest  water  surface  plane.  Water  surface  levels  for  the  Columbia  River  seldom  go 
below  this  plane.  Any  dredging  operations  for  maintaining  the  navigational 
channel  are  also  referenced  to  this  plane  and  so  any  location  on  the  river  will  have 
some  cross  sectional  area  below  this  low  water  surface.  The  Columbia  River 
Datum  can  also  be  viewed  as  the  amount  added  to  a  gauge  reading  to  reference  that 
reading  to  NGVD. 

Since  the  model  allows  for  an  offset  to  be  applied  to  any  gauge,  it  was  decided 
that  the  Astoria  gauge  values  would  be  read  into  the  model  referenced  to  MLLW. 
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Within  the  model,  an  offset  of  -3.22  feet  was  added  to  adjust  the  Astoria  gauge 
values  to  NGVD.  By  using  this  method,  NOS  data  could  be  compared  to  predicted 
tides  (also  referenced  to  MLLW)  and  any  questionable  values  could  be  easily 
located  (no  bad  values  were  found  however). 

The  observed  stage  values  acquired  are  not  as  reliable  as  the  boundary 
condition  values.  Figure  8  is  an  example  of  the  periods  with  missing  or  bad  data 
(Wauna,  Oregon  in  this  case).  Because  some  observed  stage  data  are  questionable 
at  certain  periods,  a  full  year  of  data  for  the  model  were  not  obtained  from  the 
NWRFC  office.  This  was  a  problem  for  all  of  the  upstream  stations  on  the  Lower 
Columbia  River  at  one  time  or  another.  Only  daily  high  and  low  stages  were 
gathered  for  Portland  during  the  1987  Water  Year.  Daily  high  and  low  stages 
were  gathered  for  Vancouver  prior  to  May  13, 1987.  After  that  time,  hourly  stage 
readings  were  available  for  Vancouver. 

Since  DWOPER  does  not  use  observed  upstream  station  stages  in  simulation 
computations,  questionable  observed  stage  data  will  not  affect  the  model  output. 
Because  the  computed  stages  are  independent  of  the  observed  values,  a  continuous 
series  of  computed  stages  were  used  as  comparison  values  in  the  sensitivity  analysis. 
The  set  of  comparison  data  was  computed  from  1  October  1986,  0000  hrs  to  1 
January  1987,  0000  hrs  (all  time  is  referenced  to  Pacific  Standard  Time).  These 
dates  represent  2209  simulation  hours.  Since  error-free,  continuous  records  for 
upstream  stations  was  not  possible  for  the  1987  Water  Year,  three  periods  of 
distinct  flow  conditions  were  used  in  calibrating  the  model. 
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Figure  6.  Datum  Comparisons  for  Astoria,  Oregon 


Figure  7.  Station  Relations  for  the  Columbia  River  Datum 
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Bonneville  Dam 


Figure  8.  Missing  and  Questionable  Observed  Stages  at  Wauna,  Oregon 


B.  SENSITIVITY  ANALYSIS 
1.  Friction  Term  Variations 

To  examine  how  DWOPER  output  changed  with  varying  roughness 
coefficients,  the  Manning  n  values  for  a  single  river  reach  were  increased  by  0.01 
and  a  simulation  conducted.  After  the  simulation,  the  original  Manning  n  values 
were  decreased  by  0.01  and  another  simulation  was  conducted.  In  each  simulation 
only  one  reach  had  changed  n  values,  all  other  reaches  remained  unchanged.  In 
order  to  evaluate  these  changes,  computed  values  made  with  unaltered  coefficients 
were  used  as  comparison  values.  This  was  done  to  isolate  changes  due  to  the 
different  Manning  n's  and  not  to  be  confused  with  the  inability  of  the  model  to 
match  actual  conditions.  In  all  cases,  a  constant  of  ±0.01  was  applied  to  existing 
data  (Figures  9  and  10).  In  the  case  of  decreased  n  values,  the  computed  low  water 
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stages  are  much  lower  than  the  comparison  stages  but  high  waters  changed  very 
little.  Increasing  the  roughness  coefficients  caused  low  water  stages  to  be  much 
higher  with  high  water  stages  slightly  increased.  The  summaries  presented  in 
Tables  3  and  4  show  the  results  of  these  variations  with  Figures  1 1  and  12  showing 
how  these  roughness  coefficient  changes  affect  model  results  for  Vancouver  and 
Skamokawa.  In  all  cases,  the  largest  rms  and  bias  occurred  when  the  reach  just 
downstream  of  a  station  had  altered  roughness  coefficients.  Roughness  coefficients 
changes  caused  the  greatest  modification  of  stages  just  upstream  of  the  adjusted 
reach. 

The  equations  used  to  compute  the  bias  (mean  error)  and  root  mean  square 
(rms)  are  as  follows  (Oj=  observed  stage,  Q  =  computed  stage) : 

i^-co 


The  conclusion  from  these  simulations  is  that  upstream  stations  are 
influenced  by  downstream  Manning  n  changes  but  downstream  stations  are  little 
affected  by  upstream  changes.  Any  optimizing  procedure  that  is  a  trial  and  error 
adjustment  of  the  roughness  coefficients  must  begin  at  the  most  downstream  station 
and  progress  up  the  river  system,  with  the  dynamic  tributaries  the  last  to  be 
calibrated. 
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TABLE  3.  DECREASED  MANNING  n  SIMULATIONS 
(n  values  increased  by  0.01,  values  in  feet) 


Reach 


Station 

1 

- 5 - 

5 

4 

_ 5 _ 1 

Portland 

bias 

rms 

0.110 

0.268 

0.566 

0.629 

BP9 

■TTiTTS 

mmm 

Vancouver 

0.102 

0.237 

!Hx39ii 

BsH 

HuTmH 

St  Helens 

0.1^5 

0.373 

07731 

0.817 

0.737 

0.770 

■Ml 

Longview 

0.177 

0.430 

0.769 

0.914 

-0.045 

0.227 

0.017 

0.159 

0.020 

0.094 

Wauna 

mzrm 

0.834 

1.553 

-0.073 

0.250 

BUM 

Skamokawa 

bias 

rms 

0.008 

0.037 

■H 

TABLE  4.  INCREASED  MANNING  n  SIMULATIONS 
(n  values  increased  by  0.01,  values  in  feet) 

Readi 


Station 

1 

2 

3 

4 

5 

Portland 

bias 

rms 

-Q.2ll 

0.351 

-6.707  1 
0.768 

-0.1*74  1 

0.910 

”-6.756 

0.807 

0.665 

0.092 

Vancouver 

bias 

rms 

-0.26& 

0.329 

-0.>62 

0.763 

-0.873 

0.908 

r  -0.846  ^ 

0.921 

0.6l0 

0.093 

St  Helens 

bias 

rms 

-0.364 

0.458 

-o!835  1 
0.910 

-1.012 

1.055 

-o!o57  1 

0.143 

-o!6o3 

0.060 

Longview 

bias 

rms 

-0.3^5 

0.482 

r  -0&84 
0.984 

f  0.040 
0.250 

-6.625 

0.136 

-0.010 

0.060 

Wauna 

bias 

rms 

-0.841 

1.245 

0.654 

0.179 

0.092 

-0.012 

0.048 

0.000 

0.018 

Skamokawa 

bias 

rms 

-6.661 

0.949 

ros 

0.107 

0.066 

"-6.666 

0.032 

0.001 

0.011 

0.040 


Tabular  n  +  0.01  values 


Computed  flow  (kefs) 


Figure  9.  n-Q  Relations  for  Vancouver,  Washington 


Computed  flow  (kefs) 


Figure  10.  n*Q  Relations  for  Skamokawa,  Washington 
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2.  Boundary  Condition  Variations 
a.  Astoria  tide  variations 

The  Tongue  Point  tide  gauge  near  Astoria,  Oregon  is  the  downstream 
boundary  for  DWOPER.  As  part  of  the  sensitivity  study,  the  tide  values  were 
uniformly  changed  by  ±0.5  feet  and  ±2.0  feet  Figures  13  and  14  show  how  the  ±2 
ft.  change  at  Astoria  affected  model  output  at  Vancouver  and  Skamokawa  over  a 
one  day  period.  Figure  15  summarizes  mean  values  of  the  upstream  stations  and 
how  they  are  affected  by  these  changes  in  regard  to  river  mile.  Note  that  in  Figure 
15,  the  y  axis  is  tide  difference,  not  actual  tide.  Any  change  in  the  input  tide  values 
will  influence  model  computations  far  upstream.  Model  simulations  with  the  new 
downstream  boundary  conditions  showed  that  70%  of  the  tide  difference  appears  at 
Vancouver  and  Portland,  80%  at  St.  Helens,  85%  at  Longview,  93%  at  Wauna  and 
95%  at  Skamokawa.  In  all  cases,  the  differences  are  between  previously  computed 
comparison  values  and  newly  computed  values.  These  percentages  are 
approximately  the  same  with  a  ±0.5  feet  or  ±2.0  feet  boundary  condition  change. 
This  is  due  in  part  to  the  boundary  condition  being  stage  and  not  flow  values.  Any 
increase  (or  decrease)  of  the  stage  will  considerably  change  the  amount  of  water 
entering  or  leaving  the  system.  As  the  tide  wave  travels  up  the  shallowed  sloped 
Columbia,  it  will  be  modified  by  decreased  cross  sectional  areas  and  channel 
widths.  Also  dampening  the  tide  wave  will  be  the  freshwater  runoff. 
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Figure  13.  Downstream  Boundary  Variation  Effects 
at  Vancouver,  Washington 


Time  (hrs) 


Figure  14.  Downstream  Boundary  Variation  Effects 
at  Skamokawa,  Washington 
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Figure  15.  Summary  of  Downstream  Boundary  Variation  Effects 

on  Upstream  Stations 


b.  Bonneville  Dam  discharge  variations 

The  affects  of  varying  Bonneville  flows  (the  upstream  boundary 
condition)  by  ±10%  and  ±25%  are  markedly  different  than  Astoria  tidal  stage 
changes  as  shown  in  Figures  16  and  17.  The  decreasing  bias  is  caused  by  a 
combination  of  the  increased  river  cross  sections  further  downstream  and 
proximity  to  the  downstream  boundary  condition.  Upstream,  where  the  tide 
influence  is  weakest,  the  tidal  cycle  is  more  likely  to  be  “washed  out”  by  the  higher 
flows  of  the  Columbia.  Also,  the  change  from  comparison  stages  by  these 
Bonneville  flow  variations  are  not  constant  over  time  and  fluctuate  with  the  tide 
cycle.  Figure  16  shows  a  larger  difference  for  the  -25%  Bonneville  discharge 
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variation  at  low  water  than  at  high  water  with  Figure  18  summarizing  these 
differences  for  the  entire  river.  Figure  18  also  shows  that  stage  biases  caused  by 
varying  Bonneville  outflows  are  proportional. 

As  stations  further  downstream  are  analyzed,  the  affects  of  the 
Bonneville  flow  variation  are  almost  insignificant.  In  Figure  17,  the  three 
different  curves  are  virtually  indistinguishable.  This  is  due  in  part  to  the  larger 
cross  section  areas  of  the  Columbia  River  nearer  the  mouth  and  the  proximity  of 
Skamokawa,  Washington  to  the  downstream  boundary  condition  (Astoria, 
Oregon). 


Time  (hrs) 


Figure  16.  Upstream  Boundary  Variation  Effects  at 
Vancouver,  Washington 
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Figure  17.  Upstream  Boundary  Variation  Effects  at 
Skamokawa,  Washington 


Figure  18.  Summary  of  Upstream  Boundary  Variation  Effects 
on  Downstream  Stations 


c.  Cross-section  area  variations 

The  five  cross-section  locations  closest  to  the  downstream  boundary 
were  varied  by  ±50,000  ft2  for  two  simulations  and  the  results  compared  to 
previously  computed  values  with  unchanged  areas.  As  can  be  seen  in  Figures  19 
and  20,  low  waters  are  affected  the  most  by  the  increased  areas  with  high  waters 
less  affected.  At  the  Skamokawa  gauge  site,  the  station  is  within  the  part  of  the 
river  with  changed  areas  and  so  does  not  show  the  “dam”  affect  seen  at  the 
Vancouver  gauge  site  and  other  upstream  sites.  Also  important  is  the  resultant 
decreased  amplitude  for  the  upstream  station. 


Time  (hn) 


Figure  19.  Downstream  Area  Variation  Effects 
at  Vancouver,  Washington 
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Time  (hrs) 

Figure  20.  Downstream  Area  Variation  Effects 
at  Skamokawa,  Washington 

C.  CALIBRATION  PROCEDURE 

DWOPER  has  an  automatic  calibration  feature.  This  subroutine  will  determine 
the  optimum  roughness  coefficients  which  will  minimize  the  difference  between 
observed  and  computed  stages.  This  procedure  starts  at  the  most  upstream  station 
and  progresses  downstream  reach  by  reach.  Dynamically  routed  tributaries  are 
calibrated  before  the  main  river  with  the  flows  added  to  the  main  river  as  lateral 
flows.  Upstream  and  downstream  stage  value  for  each  reach  are  needed  as  well  as 
upstream  discharge  rates.  A  steady  state  river  condition  is  assumed  by  this 
procedure. 

In  the  case  of  the  Columbia  River,  tidal  influences  make  assumptions  of  a 
steady  state  condition  invalid.  For  any  station  there  is  no  unique  stage-discharge 
relation  (as  in  the  case  of  flow  reversals).  It  is  because  of  this  complex  situation  that 
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a  trial  and  error  approach  was  taken.  The  previously  described  sensitivity  analysis 
was  made,  in  part,  to  act  as  a  guide  in  the  calibration  of  the  model. 

Simulations  for  three  periods  of  distinct  flow  conditions  were  used  in 
calibrating  the  model.  All  stations  had  error  free,  continuous  observed  records 
during  the  following  dates: 

•  21  October  1986, 1600  to  31  October  1986, 0700  (232  hours) 

•  9  January  1987, 1600  to  16  January  1987, 1600  (169  hours) 

•  8  May  1987, 1600  to  17  May  1987, 0400  (205  hours) 

The  first  period  is  a  typical  low  flow  interval  occurring  in  the  fall  with 
maximum  tide  influence  along  the  entire  course  of  the  river.  The  second  period  is 
during  a  winter  storm  runoff  event  and  is  shorter  in  duration  than  the  spring 
snowmelt  runoff  season,  represented  by  the  third  period  listed. 

In  each  case,  simulations  for  each  of  the  above  periods  were  made.  Upon 
examination  of  the  printout,  only  the  Longview  station  showed  a  peculiar  shaped 
falling  tide  hydrograph  at  different  times  during  the  May  calibration  period.  All 
other  stations  had  expected  or  typical  tide  curves  in  the  output.  Figures  21  and  22 
illustrate  the  Longview  curve  before  and  after  calibration. 

From  the  sensitivity  analysis,  most  stations  needing  calibration  were  simply  a 
matter  of  adjusting  the  appropriate  Manning  n  or  roughness  coefficient.  However, 
with  the  peculiar  situation  at  Longview,  trial  and  error  adjustment  of  the  cross 
section  area  just  above  and  just  below  the  station  was  needed.  After  Longview,  the 
river  system  was  calibrated  from  a  downstream  to  upstream  direction. 

The  adjusted  cross-sectional  areas  reflect  a  general  constriction  of  the  river  at 
this  point  although  no  data  were  available  to  verify  that  this  is  the  actual  situation. 
However,  it  was  at  this  point,  just  upstream  of  Longview  at  the  confluence  of  the 


32 


Sage  (ft)  Stage  (ft) 


Cowlitz  River,  that  debris  from  the  1980  eruption  of  Mount  Saint  Helens  entered 
the  Columbia  River. 


Time  (hn)  May  IS  and  16, 1987 


Figure  21.  Pre-Calibration  Longview 
Computed  and  Observed  Tidal  Curves 


Time  (hn)  May  IS  and  16, 1987 


Figure  22.  Post-Calibration  Longview 
Computed  and  Observed  Tidal  Curves 
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The  cross-sectional  area  variations  showed  that  the  only  way  to  increase  a 
station  high  water  value  was  to  constrict  the  river  downstream  of  the  station  of 
interest.  Varying  Manning  n  factors  would  not  affect  the  high  water  very  much  and 
boundary  conditions  were  all  screened  to  eliminate  errors  from  those  values. 
Whenever  a  station  seemed  to  need  an  adjustment  to  the  cross-sectional  areas,  the 
resulting  decrease  in  amplitude  was  increased  by  reducing  the  roughness 
coefficients  (in  effect,  smoothing  the  river  bed  to  allow  low  waters  to  be  lower).  A 
list  of  the  parameters  obtained  from  the  calibration  exercise  can  be  found  in  the 
appendix. 

The  following  tables  summarize  the  descriptive  statistics  found  for  each  of  the 
stations  for  each  of  the  different  calibration  periods  both  before  and  after 
parameter  adjustments. 


TABLE  5.  PRE-CALIBRATION  AND  POST-CALIBRATION 
BIAS  (FT)  AND  ROOT  MEAN  SQUARE  (FT) 


_ October 

Pre-  Post-  I  Pre- 
Station  calibration  calibration  calibration 


bias 

rms 


Vancouver  bias 

rms  0.483 


Portland 


Pre-  Post¬ 
calibration  calibration 


0.438 


St.  Helens 


rms  0.768 


bias 


Longview  |  rms  I  n.588 


bias 


Wauna 


Skamokawa 


rms  I  0.180 


0.205 


0.296 


0.278  0.398  0.291 


34 


D.  RESIDUAL  ANALYSIS 


The  residual  (or  stage  error)  is  the  difference  between  the  computed  (Cj)  and 
observed  stages  (0{)  as  denoted  by  Ri.  By  defining  the  residual  this  way,  any 
correction  to  the  computed  stages  will  always  be  algebraically  added  (Ri  +  Q  =  Oj). 

The  statistics  calculated  by  the  DWOPER  program  defines  the  residual  to  be 
computed  minus  observed  stage  (Q  -  Oj  =  Rj). 

Figures  23  through  28  show  a  least  squares  fit  of  observed  versus  computed 
stage  and  residual  error  versus  computed  stage.  In  no  case  were  the  correlation 
coefficients  (r2)  greater  than  0.360  for  the  residual  error  versus  computed  stage 
graphs  (most  values  between  0.00  and  0.15).  This  would  indicate  no  correlation  or 
a  non-linear  relationship  between  computed  stages  and  the  resultant  residual.  If  the 
residuals  are  random  and  normally  distributed,  it  would  be  difficult  to  improve  the 
DWOPER  output  by  modelling  these  residuals.  However,  a  time  series  of  the 
residuals  was  made  and  spectral  estimations  of  the  residuals  for  all  stations  and  all 
calibrations  are  shown  in  Figures  29  through  31 .  No  spectral  estimates  were  made 
for  Portland  or  Vancouver  because  the  observed  records  for  these  stations  were 
only  composed  of  daily  high  and  low  readings  and  not  hourly  stage  readings.  All 
other  stations  had  hourly  stage  values.  If  the  residuals  were  random  in  nature,  a 
spectral  analysis  should  show  a  white  noise  distribution  over  all  frequencies. 
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gure  25.  January  Observed  versus  Computed  Stage  Scatter  Plots 
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Figure  28.  May  Residual  versus  Computed  Stage  Scatter  Plots 
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Figure  29.  October  Residual  Spectral  Estimate 
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Figure  30.  January  Residual  Spectral  Estimate 
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FREQUENCY  FRECl/ENCY 

Figure  31.  May  Residual  Spectral  Estimate 

The  very  strong  peaks  at  specific  frequencies  indicates  that  the  residual 
components  are  strongly  tied  to  the  tidal  cycle  frequencies  of  0.081  hr1  (12.3  hrs), 
0.042  hr1  (24.0  hrs)  and  subsequent  harmonics  of  these  frequencies.  In  most  cases, 
the  12.3  hr  peak  is  dominate  and  has  the  largest  pan  of  the  residual  associated  with 

44 


1 


it.  Because  the  peak  values  vary  from  station  to  station  and  between  calibration 
periods,  the  residual  components  are  distributed  over  specific  frequencies  on  a 
seasonal  basis  for  each  site  along  the  river.  Thus,  no  one  set  of  correction  values 
can  be  applied  to  all  the  different  flow  conditions  on  the  Columbia  River. 

During  the  low  flow  October  period,  both  the  0.081  and  0.042  Hz  frequencies 
are  present  at  St.  Helens.  However,  no  one  frequency  dominates  during  the  January 
and  May  calibration  periods.  The  January  period  does  have  a  pronounced  peak 
near  the  frequency  related  with  the  length  of  the  calibration  period  itself  (169 
days).  This  may  be  associated  with  the  flow  regulation  of  Bonneville  dam  and 
Willamette  River  dams  during  the  winter  storm  event.  This  frequency  also  shows 
up  at  Longview  but  diminished  as  flow  becomes  less  a  factor  the  further  a  station  is 
downstream.  At  Skamokawa,  all  stages  for  all  calibrations  matched  extremely 
well.  The  cross  correlations  between  computed  stages  and  residuals  in  Figures  32 
to  34  show  a  sinusoidal,  not  linear,  correlation  exists. 

A  one  to  three  hour  phase  shift  between  the  observed  stages  and  the  stage  error 
is  present  in  each  of  the  calibration  periods  for  all  stations.  The  phase  shift  becomes 
evident  when  an  observed  stage  time  series  is  back-lagged  with  respect  to  the 
corresponding  stage  error  time  series.  All  the  stations,  except  for  Skamokawa, 
have  similarly  shaped  cross  correlation  curves  with  respect  to  each  other,  especially 
during  the  October  calibration  period.  In  each  of  the  calibration  periods, 
Skamokawa  is  always  180°  out  of  phase  with  the  other  stations,  except  for 
Longview  during  the  May  calibration  period.  All  stations  had  larger  correlations 
(either  negative  or  positive)  at  the  12  to  13  hour  and  24  to  25  hour  forward  and 
back-lag  times. 


It  should  be  noted  that  for  the  October  calibration,  the  observed  and  computed 
lower  low  water  times  for  all  stations  matched  extremely  well  (0  to  1  hour 
difference).  For  Skamokawa,  times  for  computed  and  observed  high  and  low 
waters  nearly  always  matched  no  matter  which  calibration  period  was  examined.  It 
is  evident  that  even  though  DWOPER  handles  the  timing  of  the  extreme  parts  of  the 
tide  cycle  well,  intermediate  stages  are  not  modeled  as  well. 
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Figure  34.  Cross  Correlations  for  May 

The  Northwest  River  Forecast  Center  is  presently  simulating  river  stages  for 
the  three  days  prior  to  a  forecast  period  (a  three  day  “backup  period”)  as  well  as  the 
forecast  period  itself.  An  analysis  of  how  well  the  computed  stages  are  reproducing 
actual  stages  is  possible  by  having  a  backup  period.  The  backup  period  data  also 
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allows  for  a  method  of  adjusting  the  computed  stages  based  on  the  most  recent  river 
conditions.  Adjustment  factors  are  the  average  differences  between  the  observed 
and  computed  stages  for  four  stage  ranges.  The  adjustment  factors  are  trended  into 
the  future  and  are  used  as  corrections  to  the  calculated  river  stages.  Each  stations  is 
corrected  separately  with  the  forecasted  stages  placed  in  the  CROHMS  regional 
database  for  outside  users  (Orwig,  1986). 

It  has  been  shown  that  the  residuals  occur  at  the  same  frequency  as  the  tide 
cycle  and  implies  that  as  the  tide  changes,  so  will  the  residuals.  Figure  35  illustrates 
how  the  residuals  vary  with  the  computed  stage  over  a  single  tide  cycle.  The 
Longview  looped  error  curve  was  developed  from  the  October  calibration  period. 
The  first  part  of  the  curve  (points  1  to  2)  is  from  the  higher  high  to  lower  low  tide, 
the  second  part  (2  to  3)  is  from  lower  low  to  lower  high,  the  third  part  (3  to  4)  is 
from  lower  high  to  higher  low  and  the  fourth  part  (4  to  5)  is  from  higher  low  to  the 
higher  high  stage  of  the  next  tide  cycle. 

For  any  one  stage  value,  the  adjustment  factor  to  be  added  to  the  computed 
stage  will  depend  on  what  part  of  the  tide  cycle  that  stage  falls.  The  current 
procedure  at  NWRFC  does  not  take  this  tide  cycle  variation  into  account  but  such  a 
procedure  may  prove  useful  at  higher  river  flow  periods  where  tidal  influences  are 
weakened  or  washed  out  altogether. 

The  looped  error  curves  vary  in  shape  from  station  to  station  but  are  similar  at 
a  station  for  comparable  tide  cycles.  As  the  tide  cycle  gradually  changes  (when  the 
lower  high  becomes  the  higher  high  for  example)  the  error  curve  changes  shape  as 
well.  Any  adjustment  scheme  using  looped  error  curves  will  still  need  some  sort  of 
a  backup  period  analysis  to  allow  error  curves  to  be  based  on  the  most  recent  river 
conditions.  An  historic  record  of  curves  for  differing  tidal  cycles  at  a  station  can  be 
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be  obtained  as  these  curves  are  generated.  When  forecasted  river  flow  conditions 
change  markedly  from  the  backup  period  or  if  observed  data  is  not  available,  any 
backup  period  curves  become  useless  and  the  historic  record  of  curves  could  be 
used  to  find  appropriate  adjustment  factors. 

In  the  case  of  Portland  and  Vancouver,  the  method  described  above  cannot  be 
applied  because  only  high  and  low  tide  stages  were  recorded.  To  have  an  accurate 
and  functional  adjustment  procedure,  hourly  stage  data  must  be  available.  Any 
plan  to  decrease  the  amount  of  data  gathered  at  any  of  the  lower  Columbia  River 
sites  would  lessen  a  hydrologist's  ability  to  produce  an  accurate,  knowledgeable 
forecast  for  the  river. 
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IV  MODEL  APPLICATIONS 


The  main  use  of  the  DWOPER  program  has  been  to  forecast  stages  on  the 
Lower  Columbia  River  to  support  the  shipping  interests  in  the  area.  As  outlined 
earlier  in  this  paper,  the  Port  of  Portland  and  other  port  authorities  on  the  lower 
Columbia  River  have  used  these  forecasts  to  increase  the  amount  of  commerce  on 
the  river.  However,  other  possible  uses  of  the  DWOPER  model  are  discussed 
below. 

A.  FLOOD  ROUTING 

A  primary  reason  the  National  Weather  Service  produced  DWOPER  was  for 
the  routing  of  streamflow  on  major  river  systems.  As  such,  this  model  has  been 
used  by  the  U.S.  Army  Corps  of  Engineers  to  simulate  downstream  river  stages 
with  various  Bonneville  dam  regulation  schedules.  By  doing  so,  downstream 
effects  can  be  approximated  and  potential  flood  damage  avoided.  On  the  other 
hand,  estimated  dam  regulation  for  maintaining  minimum  river  flow  conditions 
during  the  dry  summer  and  autumn  periods  has  been  simulated  by  using  the 
DWOPER  program. 

Also  important  is  the  link  the  Northwest  River  Forecast  Center  (NWRFC)  has 
to  the  Corps  of  Engineers  Columbia  River  Operational  Hydromet  Management 
System  (CROHMS)  regional  database.  By  using  the  Corps'  Synthetic  Streamflow 
and  Reservoir  Regulation  (SSARR)  watershed  and  river  routing  model  as  input  to 
the  DWOPER  program,  the  latest  and  best  estimates  of  lateral  inflows  and  ocean 
tides  can  be  used  to  produce  simulations  on  an  operational  basis  (Orwig  et  al., 
1986). 
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B.  CONFIDENCE  INTERVALS 


The  NWRFC  has  successfully  used  the  DWOPER  model  in  producing  three  day 
Lower  Columbia  stage  forecasts  since  May,  1986.  This  office  has  now  been  asked 
by  the  local  merchants  exchange  to  extend  the  river  stage  forecasts  from  three  days 
to  six  days  (personal  conversation  with  Orwig).  Questions  of  input  data  quality 
control  and  output  adjustments  must  be  addressed  when  forecasting  stages  so  far 
into  the  future.  Errors  are  a  function  of  model  exactness  and  input  data  accuracy. 

The  sensitivity  analysis  of  DWOPER  described  earlier  in  this  study  can  be  used 
by  forecasters  to  quantify  how  boundary  condition  variations  will  affect  model 
output.  By  plotting  computed  stages  made  under  varying  Bonneville  discharges,  a 
type  of  “confidence  interval”  for  any  station  can  be  produced.  Such  an  approach 
would  help  estimate  the  possible  river  stages  cause  by  unexpected  changes  in 
Bonneville  dam  releases. 

In  Figure  36,  the  effects  of  changing  Bonneville  dam  outflows  by  ±25%  at 
Portland  are  shown.  The  low  river  flow  during  this  time  had  a  large  tidal 
influence.  The  lower  estimated  flow  (-25%)  produced  lower  stages  overall  with 
the  low  tide  having  a  slightly  larger  offset  than  the  high  water  tide  from  the 
comparison  stages.  The  higher  estimated  flow  (+25%)  produced  higher  overall 
stages  with  low  water  depths  being  modified  the  most.  The  differences  over  this 
tide  cycle  are  depicted  in  Figure  37.  This  graph  indicates  that  as  Bonneville  dam 
discharges  are  overestimated,  forecasts  will  have  a  greater  error  at  the  low  waters. 
When  dam  discharges  are  underestimated,  an  offset  of  approximately  0.6  feet  is 
produced  for  all  stages. 
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Figure  36.  Bonneville  Discharge  Variation  Effects 
at  Portland,  Oregon  (October  5-6,  1986) 


Figure  37.  Stages  Differences  From  Bonneville  Discharge  Variation  at 
Portland,  Oregon  (October  5-6,  1986) 


In  Figure  38,  the  effects  of  changing  Bonneville  dam  outflows  by  ±25%  at 
Portland  during  a  high  river  flow  period  is  shown.  During  this  time,  there  was  less 
tidal  influence  than  in  the  October  period.  The  lower  estimated  flow  (-25%)  again 
produced  lower  stages  overall.  The  higher  estimated  flow  (+25%)  produced 
higher  overall  stages  with  the  second  low  water  almost  eliminated  (Figure  39).  The 
graph  indicates  that  as  Bonneville  dam  discharges  are  overestimated,  forecasted 
stages  become  very  uniform  with  tidal  influences  almost  entirely  removed.  The 
result  will  be  greater  errors  occurring  at  the  low  waters.  When  dam  discharges  are 
underestimated,  a  false  tidal  influence  could  be  introduced  where  none  exists.  In 
the  high  river  flow  case,  departures  from  the  comparison  stages  are  much  less  tied 
to  the  tide  cycle  and  are  more  complex. 


Figure  38.  Bonneville  Discharge  Variation  Effects 
at  Portland,  Oregon  (November  25,  1986) 
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Figure  39.  Stages  Differences  From  Bonneville  Discharge  Variation  at 
Portland,  Oregon  (November  25,  1986) 


The  concept  of  over  and  underestimated  boundary  conditions  can  be  applied  to 
the  Tongue  Point  tide  gauge  as  well.  When  a  forecaster  estimates  a  storm  surge 
influence  on  the  ocean  tides,  the  forecaster  will  know  how  far  upstream  the  surge 
will  be  felt.  Figures  39  and  41  show  how  Skamokawa  stages  change  by  adjusting 
the  Astoria  tide  values  by  ±  2  feet.  During  this  period,  tidal  influence  on  the 
Skamokawa  stages  is  very  evident.  In  both  periods,  the  higher  estimated  tide  (+  2 
ft)  produced  higher  overall  stages  while  the  lower  estimates  produced  lower 
overall  stages.  Differences  over  the  tide  cycles  are  depicted  in  Figures  40  and  42. 
Figure  41  indicates  that  the  departures  from  the  comparison  stages  are  not  just  a 
constant  offset  and  that  the  departures  are  associated  with  the  tide  cycle. 
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Figure  40.  Astoria  Tide  Variation  Effects  at  Skamokawa,  Washington 

(October  5-6,  1986) 


Figure  41.  Stages  Differences  From  Astoria  Tide  Variation 
at  Skamokawa,  Washington  (October  5-6,  1986) 
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Figure  42.  Astoria  Tide  Variation  Effects  at  Skamokawa,  Washington 

(November  25-26,  1986) 


Figure  43.  Stages  Differences  From  Astoria  Tide  Variation 
at  Skamokawa,  Washington  (November  25-26,  1986) 
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C.  OTHER  MAJOR  USES 

Normally  the  DWOPER  program  provides  output  such  as  graphs  displaying 
computed  and  observed  values  at  gauging  stations,  a  list  of  numbers  plotted  and  the 
statistics  of  bias  and  root  mean  square.  Additionally,  the  following  three  special 
tabular  printouts  can  be  specified. 

•  Computed  stages,  discharges  and  velocities,  and  observed  stages  for  all 

stations  on  one  river  can  be  displayed  for  one  time  step. 

•  Computed  and  observed  stages  at  a  single  station  can  be  displayed  for  a 

period  up  to  25  hours. 

•  Computed  discharges  and  velocities  for  a  single  station  can  be  displayed 

for  a  period  up  to  25  hours. 

The  ability  to  compute  average  cross  sectional  velocities  can  be  helpful  in 
predicting  average  river  currents.  Applications  could  include  modelling  oil  spill 
spreading,  special  boat  operations,  fish  migration  or  location  of  objects  lost  in  the 
river. 

Tidal  range  variations  influenced  by  dredging  can  be  studied  by  varying  cross 
sectional  areas.  By  decreasing  the  cross  sectional  areas,  debris  or  obstructions 
influences  in  the  river  can  be  studied. 

The  DWOPER  model  can  be  used  as  an  error  check  since  the  program  does  not 
use  observed  stages  in  computing  stages  for  a  station.  With  the  calibration  of  the 
model,  residuals  were  kept  to  within  1.5  feet  of  the  observed  stages.  By  using  the 
model  output  as  a  estimator  for  observed  stages,  stage  data  could  be  screened  before 
entry  into  the  Columbia  River  Operational  Hydromet  Management  System 
(CROHMS)  regional  database.  This  database  is  a  key  link  between  several  federal 
agencies  and  other  public  and  private  users  in  the  Pacific  Northwest  (Pasteris, 
1984). 
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V.  SUMMARY  AND  CONCLUSIONS 


The  DWOPER  model  uses  a  weighted  four-point  implicit  finite  difference 
method  to  solve  the  Saint  Venant  flow  equations  of  conservation  of  mass  and 
momentum.  Both  upstream  and  downstream  boundary  conditions  of  river  flow  or 
stage  are  specified  to  obtain  solutions  to  h  and  Q  along  intermediate  points  on  the 
river.  The  Dynamic  Wave  model  is  a  mathematical  model  well  suited  for  a  variety 
of  rivers  and  flow  conditions.  As  applied  to  the  Columbia  River,  this  model  can 
provide  essential  information  on  hourly  river  stages  critical  for  navigational  safety 
and  ship  commerce  on  this  economically  important  river. 

A  sensitivity  analysis  was  performed  by  varying  the  roughness,  boundary 
conditions,  and  cross-sectional  areas.  The  computed  values  from  1  October  1986, 
0000  hours  to  1  January  1987,  0000  hours  were  used  as  comparisons.  These  dates 
represent  2209  simulation  hours.  The  calibration  procedure  examined  the  distinct 
flow  conditions  of  a  low  river  flow  period,  a  winter  storm  event,  and  a  spring 
snowmelt  season. 

In  testing  the  Manning  n  roughness  coefficients,  each  set  of  Manning  n  values 
were  varied  by  ±0.01  for  each  of  the  river  reaches.  In  the  case  of  decreased  n 
values,  the  computed  low  water  stages  were  much  lower  than  the  comparison 
stages;  high  waters  changed  very  little.  Increasing  the  roughness  coefficients 
caused  low  water  stages  to  be  much  higher  with  high  water  stages  slightly 
increased.  The  conclusion  from  these  simulations  is  that  upstream  stations  are 
influenced  by  downstream  Manning  n  changes  but  downstream  stations  are  little 
affected  by  upstream  changes.  Any  optimization  procedure  that  is  a  trial  and  error 
adjustment  of  the  roughness  coefficients  must  begin  at  the  most  downstream  station 
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and  progress  up  the  river  system,  with  the  dynamic  tributaries  the  last  to  be 
calibrated. 

Downstream  boundary  changes  of  ±0.5  feet  and  ±2.0  feet  were  made  to  the 
Astoria  tide  stages.  Model  simulations  showed  that  70%  of  the  tide  difference 
appears  at  Vancouver  and  Portland,  80%  at  St.  Helens,  85%  at  Longview,  93%  at 
Wauna  and  95%  at  Skamokawa.  These  differences  occurred  as  a  constant  offset  and 
did  not  appear  to  be  associated  with  the  tide  cycle. 

The  effects  of  varying  the  upstream  boundary  condition  (Bonneville  Dam 
discharges)  by  ±10%  and  ±25%  were  markedly  different  from  the  downstream 
boundary  changes.  Upstream,  where  the  tide  influence  is  weakest,  the  tidal  cycle  is 
more  likely  to  be  “washed  out”  by  the  higher  flows  of  the  Columbia.  Also  these 
changes  fluctuated  with  the  tide  cycle.  Downstream  stations  did  not  show  such 
differences  because  of  the  larger  cross  section  areas  of  the  Columbia  River  nearer 
the  mouth  and  the  proximity  to  the  downstream  boundary  condition. 

To  observe  the  effects  of  varying  cross-sectional  areas,  the  five  most 
downstream  cross-section  locations  were  changed  by  ±  50,000  ft2.  Low  waters 
were  changed  the  most  by  the  increased  areas  with  high  waters  less  influenced.  At 
the  Skamokawa  gauge  site,  the  station  is  within  the  part  of  the  river  with  altered 
areas  and  so  did  not  show  the  “dam”  affect  seen  at  the  Vancouver  gauge  site  and 
other  upstream  sites.  Because  the  tide  wave  could  not  easily  propagate  upstream 
with  the  restricted  cross  section  areas,  a  decreased  tide  amplitude  for  upstream 
stations  occurred.  Fresh  water  runoff  encountered  the  decreased  areas  as  a  “dam” 
and  backed  up  causing  an  overall  increase  in  upstream  stages. 

Most  stations  needing  calibration  involved  a  simple  matter  of  adjusting  the 
appropriate  Manning  n  or  roughness  coefficient.  The  adjusted  cross-sectional  areas 
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areas  reflect  a  general  constriction  of  the  river  at  this  point,  although  no  data  were 
available  to  verify  that  this  is  the  actual  situation.  However,  it  was  at  this  point,  just 
upstream  of  Longview  at  the  confluence  of  the  Cowlitz  River,  that  debris  from  the 
1980  eruption  of  Mount  Saint  Helens  entered  the  Columbia  River.  As  part  of  the 
Longview  calibration  two  extra  cross  sections,  one  just  above  and  another  just 
below  the  Cowlitz  River  confluence,  were  added.  After  the  Longview  adjustments, 
the  river  system  was  calibrated  from  a  downstream  to  upstream  direction. 
Comparisons  of  observed  data  to  post-calibration  simulations  showed  that  for  the 
October  and  January  periods,  five  of  six  stations  had  consistently  lower  mean  stage 
errors  while  in  the  May  period,  three  of  six  stations  had  consistently  lower  mean 
stage  errors. 

Regression  analyses  of  the  computed  residual  values  for  each  of  the  stations 
gave  correlation  coefficients  (r2)  less  than  0.360.  However,  cross  correlations 
between  residual  and  computed  stages  showed  that  the  two  series  were  highly 
sinusoidally  correlated  for  all  stations.  A  spectral  estimation  of  the  residuals 
exhibited  very  strong  peaks  at  frequencies  of  0.081  hr1  (12.3  hrs),  0.042  hr1  (24.0 
hrs)  and  subsequent  harmonics  of  these  frequencies.  The  residual  components  are 
strongly  associated  with  the  tidal  cycle. 

The  spectral  values  of  the  residuals  changed  between  calibration  periods  and 
indicates  that  river  flow  condition  influences  the  residual  distribution.  A  one  to 
three  hour  phase  shift  between  the  observed  stages  and  the  stage  error  is  present  in 
each  of  the  calibration  periods  for  all  stations.  All  stations  had  larger  correlations 
(either  negative  or  positive)  at  the  12  to  13  hour  and  24  to  25  hour  forward  and 
back-lag  times  from  the  cross  correlation  peak. 
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The  DWOPER  program  can  be  used  in  the  areas  of  forecasting  stages, 
estimating  average  velocities,  effects  of  dredging  or  constrictions  due  to  added 
debris.  As  this  type  of  information  is  used  by  the  public,  greater  demands  for 
longer,  more  accurate  forecasts  will  be  made.  The  NWRFC  project  to  extend 
forecasts  to  six  days  is  such  an  example. 

The  study  presented  here  can  be  used  for  a  basis  in  understanding  how  the 
DWOPER  program  reacts  to  changes  in  boundary  conditions  and  how  to  approach 
predicting  stage  errors.  Necessary  adjustment  of  model  results  must  be  made  to 
provide  even  higher  quality,  useful  data.  Ideas  like  confidence  intervals  for  stages 
based  on  estimations  of  dam  discharge  variability  can  be  addressed  by  simulations 
made  beforehand. 

Also  demonstrated  by  this  study  is  the  need  for  dependable,  high  quality  stage 
measurements.  Any  type  of  short  term  corrections  to  model  output  will  need 
recent  accurate  hourly  records  for  all  stations  of  interest.  Accurate  hourly  records 
will  allow  a  correction  algorithm  to  be  used  based  on  the  dominant  frequencies.  As 
these  frequencies  change  with  flow  conditions,  so  must  the  analysis  of  the  residuals. 
This  study  looked  briefly  at  three  period  of  about  200  hours  (606  hours  out  of  a 
possible  8760  hours  for  a  year)  and  can  provide  only  some  historical  reference 
data.  Further  research  is  needed  to  completely  understand  how  residuals  change 
over  time  and  over  transitional  flow  conditions.  The  1987  Water  Year  did  provide 
an  excellent  chance  to  study  DWOPER  over  a  period  when  low  stages  are  critical. 

DWOPER  handles  the  timing  of  the  extreme  parts  of  the  tide  cycle  well  but 
intermediate  stages  are  not  being  modeled  as  well.  Further  studies  of  these  residuals 
is  needed  to  develop  adjustment  procedures  to  improve  DWOPER  model  forecasts. 


APPENDIX 


DWOPER  PARAMETER  LISTING 
The  following  DWOPER  program  output  is  a  listing  of  the  carryover  file 
parameters.  These  values  describe  such  things  as  bed  roughness  coefficients  and 
cross  section  area  locations.  For  a  detailed  explanation  of  each  of  the  variables,  the 
reader  should  refer  to  the  model  documentation  (Fread,  1978). 
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